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Abstract. A large system of ordinary differential equations is approximated 
by a parabolic partial differential equation with dynamic boundary condition. 
Using the theory of differential operators with Wentzell boundary conditions, 
we give estimates on the order of approximation. The theory is demonstrated 
on a voter model where the Fourier method applied to the PDE seems to be 
of great advantage. 



1. Introduction 

It has been known for a long time that there is a wonderful interplay between 
discrete time stochastic processes and partial differential equations, see the seminal 
paper by Courant, Friedrichs and Lewy [BJ. Since then, the pioneering work of 
Feller revealed the deep connections to second order differential equations with 
"complicated" boundary conditions, see the monograph by Mandl [TT] . 

Our intention with this work is to go back to the roots and explore the con- 
nections of large systems of ordinary differential equations to parabolic partial dif- 
ferential equations with so-called "Wentzell" boundary conditions from a rather 
particular point of view: given a large system of ordinary differential equations, we 
construct an "approximating" partial differential equation, give estimates on the 
accuracy of this approximation, and show that in some cases it is much easier to 
handle the parabolic equation than the large ODE system. 

The main motivation of this theoretical investigation is to approximate a dy- 
namic process on a network by a partial differential equation. The network is given 
by an undirected graph and the process is specified by the possible states of the 
nodes and the transition probabilities. Typical examples are epidemic processes 
and opinion propagation on networks. Analysing the mean field approximation for 
the expected number of infected nodes in an epidemic process on a large network 
we were led to a first order PDE approximation in our previous work Batkai et al. 
[5]. In a recent paper, in which regular random, Erdos-Renyi, bimodal random, and 
Barabasi- Albert graphs are studied, we have shown that a suitable choice of the 
coefficients in the master equation leads to an ODE approximation with tridiagonal 
transition rate matrices (see Nagy, Kiss and Simon [I!])- Our study in this paper 
aims at approximating dynamic processes on networks, for which the transition rate 
matrix of the underlying Markov chain has a tridiagonal structure (with a possible 
extension to similar matrices). 

The paper is organized as follows. First, we introduce our general notation and 
setup along with a standard heuristic derivation of the approximating PDE via 
finite differences. Then we present the operator semigroup theoretic setup and 
the needed general approximation theorems. In Section 4 we show how to use the 
well-developed operator matrix approach to differential operators with Wentzell 
boundary conditions due to Engel and coauthors [3 |51 0] to this problem and 
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prove error estimates. Finally, in the last section we illustrate our results with two 
examples: The first one is the propagation of two opinions along a cycle graph, 
called a voter-like model, the second is an SIS type epidemic propagation on a 
complete graph. 



2. Setup 

In this section we fix our notation, collect the main definitions, derive the approxi- 
mating partial differential equations, and give the main heuristics which lie behind 
our approximation. 

Let N £ N be a large, fixed integer, and a, b and c real- valued functions on 
[-jf, 1 + jj] For < k < N, let a k := a(-|), b k := &(-£) and c k := c(-|). Consider 
the following tridiagonal matrix 
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and the corresponding (ODE) system 

f x{t) = A N x{t) 



(1) 



j(0) = ve C N+1 



on C" +1 . 

We wish to approximate the solution x(t) to this (ODE) by considering it as a 
discretisation of a continuous function u(t, z) on the interval [0, 1], i.e. 

k 



u t, 



N 



x k (t) 



for < k < N. We now derive an approximate (PDE) for the function u(-, ■) using 
the (ODE) given above. For any 1 < k < N — I we have: 

k 



d t u ( t, — 



= x k (t) = a fc _ia; fc _i(t) + b k x k {t) + c k+1 x k+l {t) 
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By considering the approximations 
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using the functions a, 6 and c, and writing h := we obtain the approximate 
(PDE) 

ft 2 



(2) 



9 t u(i, z) ~ — (a (z - h) + c (z + ft)) d zz u(t, z) 
+ h(c(z + h) — a(z — h))d z u(t, z) 
+ (a(z -h) + b(z) + c(z + h))u(t, z), 



valid for z £ (0, 1). Note that the approximation is of order 

On the boundary, similar transformations yield the first order boundary equa- 
tions 

(3) d t u(t, 0) ~ hc(h)d z u{t, 0) + (c(h) + b(0))u(t, 0) 
and 

(4) d t u{t, 1) ~ -ha(l - h)d z u{t, 1) + (<z(l - h) + b(l))u(t, 1). 

Note that here the approximations are only of order -A^. 

The initial condition u(0, z) is to be chosen as a suitable interpolation of the 
values x fc (0) at z = ^ (0 < k < N). 

In the following we consider the exact PDE and its solution, the latter being 
the approximation of the exact solution to the ODE at the points -j™, and show 
estimates on how good this approximation is. 



3. THEOREMS 

Now we give a rather general setup to prove the desired estimates on the approx- 
imation. We use the theory of operator semigroups and our general reference is 
Engel and Nagel [5]. 

Assumptions 3.1. Let X n , X be Banach spaces and assume that there are 
bounded linear operators P n : X — > X n , J n : X n — > X with the following properties: 

• There is a constant K > with ||P n ||, ||J„|| < K for all n £ N, 

• P n J n = I n , the identity operator on X n , and 

• JnPnf —> f as n — > oo for all f £ X. 

Assumptions 3.2. Suppose that the operators A n , A generate strongly continuous 
semigroups on X n and X, respectively, and that there are constants M > 0, u> £ K 
such that the stability condition 

(5) ll^n(*)ll < Me ut holds for all n £ N, t > 0. 

We will make use of a special variant of the Trotter-Kato theorem, which we cite 
here for convenience, see the lectures by Batkai, Csomos, Farkas and Ostermann 
2, Proposition 3.8]. 
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Proposition 3.3. Suppose that Assumvtions \3. I\ hold, that there is a dense subset 
Y C D (A) invariant under the semigroup T such that P n Y C D(A n ), and that Y 
is a Banach space with some norm \\ ■ \\y satisfying 

\\T{t)\\ Y < Me"*. 

// there are constants C > and p£ N with the property that for all f G Y 

\\A n P n f - P n Af\\ Xn < C^|^, 
then for each t > there is C' > such that 

\\T n (t)P n f - P n T(t)f\\ Xn < C'tt. 
Moreover, this convergence is uniform in t on compact intervals. 

This result can be slightly improved in case analytic semigroups are involved. 



Lemma 3.4. Suppose that the conditions of Proposition \373\ are satisfied and that 
A generates an analytic semigroup. If there is e £ (0, 1) and there are spaces 
Y ^ D(A) ^ Z ^ X such that T{s)Z C Y for all s > and 

\\T(s)\\ c(Z ,Y) < 

holds, then 

vll/IU 



\\T N (t)P N f-P N T(t)f\\ XN <C 



NP 



Note that this condition is for example satisfied if there is a € (0, 1) so that 
Y = D((I - ^) 1+a ) and D(A) ^ Z D((I - A) a+e ) holds. 

Proof. As in the proof of Batkai, Kiss, Sikolya and Simon (5j Lemma 5], we have 
the representation 

{P N T{t)-T N {t)P N )f = f T N {t-s)(P N A-A N P N )T{s)fds 

Jo 

for all / € D(A). Hence, using the analyticy of the semigroup T, we obtain the 
norm estimate 

\\P N T(t)f-T N (t)P N f\\< J Me^ t - s, >\\(P N A-A N P N )T(s)f\\ds 



Np Np 



where the constants M', C and M" only depend on t > 0. □ 

We have seen in the calculations of the previous section that our approximation 
is of third order in the interior of the interval and of second order on the boundary. 
Let us formalize now the calculations and put them into the general framework 
presented above. 

Theorem 3.5. Consider the ordinary differential eguation given by (TTJ) and the 
approximating partial differential equation ([2]) with dynamic boundary conditions 
Q and (gl), where v = Pnu(0, •)• If there is e € (0, \) such that u(0, •) e D((I - 
A)i +£ ), then for all T > there is C = C{T) > such that for all t G (0, T] we get 

(6) \\P N u(t, •) - as^Hoo < ^3 ||u(0, ■)\\ D{{I _ A)i+r 
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Proof. As a first step, we have to associate to the partial differential equation 
@ with boundary conditions © and (HJ) a Banach space X and a generator A. 
Following the treatment in Engel [7j [8] or in Batkai and Engel [4] , we introduce the 
spaces 

X:=C[0,1], 

and 

X :={(/) eX xC 2 |y=(.f(0),/(l)) T }. 
Let us also consider the operators 

(D m f)(z) := y (a (* - h) + c(z + h)) f"(z) 

+ h{c{z + h)- a(z - h))f{z) + (a(z - h) + b(z) + c(z + h))f(z) 
defined on its maximal possible domain, D(D m ) := C 2 [0, 1], and 

Uf. = ( hc(h)f'(0) + (c(h)+b(0))f(0) \ 
3 ' \-ha(l-h)f'(l) + (a(l-h)+b(l))f(l) J 

defined on D(D m ) and mapping to C 2 . Our associated operator should be 



Af = D m f with D{A) := {/ G D(D m ) : ( D m f(o),D m f(i) ) T = Bf} . 

Then, by Engel [7] or Batkai and Engel [H Remark 4.4], this operator generates an 
analytic semigroup of angle ^ in the space X. 
Further, we introduce the spaces 

: JV - 1 x c 2 



and define the operators P, 



N 



Xn '■= 

: X — > Xn as 



P N (f,y)~{x,y) with x k = f(±) forfe = l,2,...,JV-l, 
and take J/v to be a suitable interpolation. Clearly, these operators and spaces 
satisfy the conditions in Assumptions 13.11 

Abusing the matrix notation, define now on Xn — C N+1 the operator 
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Taking (J) G £>(A) (i.e., / € W 1 ' 1 ^, 1) and yi = /(0), y 2 = /(l), we see that 

00/(0) + fc/OO+ca/^) 

ai/(&)+&a/($) +<*/(&) 



'A - 



aN-2j(- 
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6o/(0) + Cl /(i) 
a N _ 1 /(^i) + 6^/(1) 
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and that 



+ A7 W^m^AT) - c(^)/(^))/'(#) 

+ W^) + K#) + c(^)/(£) 
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for k = 1,2,..., N- 1, and 

(P N A(f)) N =±c(±)f(0) + (c(i) + 6(0))/(0) 
(P N A(f)) N+1 = - i c (^i)/'(l) + (a(^) + 6(1))/(1) 
By the calculations of the previous section, we see that there is C > such that 



\(p N A{f y )) k - (A N p N (f))k\ <j^\\r\u 

\(P N A({)) N - (A N P N {l))N\ <^||/"||oo, 
\(P N A(f)) N+ i - (A N P N (f y )) N+1 \ <^||/"||oo 

holds. Since A generates an analytic semigroup, it leaves Y = C 3 [0, 1] invariant. 
Hence, Proposition 13.31 is applicable with Y — C 3 [0, 1] and we obtain the desired 
estimate for all u(0, •) £ Y. To improve on this result and relax the regularity 
assumption on the initial value, we use the analyticity of the semigroup and Lemma 
13 with a = \. 

Introducing the notation B = I — A, our aim now is to show that D(B 3 / 2 ) c 
C 2 [0,1] n C 3 (0,1). Since D(B) C C 2 [0,1], it is enough to show that D{B^ 2 ) C 
C X (0,1). Let / G D{B 1 / 2 ) such that g = B 1 / 2 /. Then by Engel and Nagel [9, 
Corollary II.5.28], 



f°° 1 

/ -=R(X+l,A)gdX. 
Jo v X 



Further, by checking the explicit representation of the resolvent as in the proof 
in Engel and Nagel [5J Theorem VI. 4. 5], we see that the resolvent is given by the 
combination of exponential terms and a convolution term. Since we are in the 
interior of the domain, we can drop the exponential terms because they do not 
disturb regularity and concentrate on the convolution term. Hence we may assume 
that 

,oo i ^ f 1 e-^+i\-*\g( s )dsdX. 

! VA2 % /ATUo 



ii 



Rewriting, we obtain 



/(*) = / - \ / e-^^gis) ds 
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* ( e -VA+T- f e^sgis) ds 
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Formally differentiating with respect to x behind the first integral, we obtain 
' 1 " 1 -.e~^ X + Ix f e^ X + Is g{s)ds + g(x) 



o 2V^VATT I VX + 1 
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Since this improper integral converges uniformly in x on any closed subinterval 
of (0,1), and depends continuously on x, the function / is indeed continuously 
differentiable on (0, 1). 

□ 



4. Applications 

The main motivation of the previous theoretical investigation is to approximate 
a dynamic process on a network with a partial differential equation and to justify 
empirical observations. The network is usually given by an undirected graph and 
the process can be specified by the possible states of the nodes and the transition 
rate probabilities. The latter is the probability that the state of a node changes 
from one state to another depending on the states of the neighbouring nodes. In 
certain classes of models the complete state space can be reduced (using e.g. mean 
field approximations or structural symmetries), leading to tridiagonal systems. 

In this section we show in two cases how the theory can be applied. The first one 
is the propagation of two opinions along a cycle graph, called a voter-like model, 
the second is an SIS type epidemic propagation on a complete graph. As usual 
in the literature, in both models the natural Markov process is conditioned on not 
reaching the absorbing state(s). 

4.1. Voter-like model on a cycle graph. Let us consider a cycle graph with 
N +3 nodes, i.e. we have a connected graph, in which each node has two neighbours. 
A node can be in one of two states, let us denote them by and 1. These states 
represent two opinions propagating along the edges of the graph (see Holley and 
Liggett [IH]). If a node is in state and has k neighbours in state 1 (k = 0, 1, 2), 
then its state will change to 1 with probability krAt in a small time interval At. 
This expresses that opinion 1 invades that node. The opposite case can also happen, 
that is a node in state 1 can become a node with opinion with a probability kjAt 
in a small time interval At, if it has k neighbours in state 0. The parameters r and 
7 characterize the strengths of the two opinions. The model originates in physics, 
where in a network of interacting particles each node holds either spin 1 or -1 (see 
Vazquez and Eguiluz [13] )■ In a single event, a randomly chosen node adopts the 
spin of one of its neighbors, also chosen at random. 

Assuming that at the initial instant the territories of the two opinions are con- 
nected sets, the underlying conditioned Markov chain can be given as follows. The 
state space is the set {0, 1, 2, . . . , N}, where a number k represents the state in 
which there are k + 1 nodes in state 1 and they form a connected arc along the 
cycle graph. Starting from state k the system can move either to state k + 1 or 
to k — 1, since at a given instant only one node can change its state (by using 
the usual assumption that the changes at the nodes can be given by independent 
Poisson processes). When the system moves from state k to k + 1 then a new node 
in state 1 appears at one of the two ends of the arc of state 1 nodes. Hence the 
rate of this transition is 2r, expressing that a node in state and having a single 
neighbour in state 1 becomes a state 1 node, and this can happen at both ends of 
the state 1 territory. Similarly, the rate of transition from state k to k — 1 is 2"f. 
Let us denote by Xk(t) the probability that the system is in state k. The above 
transition rates lead to the differential equation 



x(t) = 2rx k -x{t) - 2(r + j)x(t) + 2jx k+1 (t). 
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(For k — and for k = N the equations contain only two terms.) Thus our system of 
ODEs takes the form given in ([T]) with a = 2r, c = 27 and 6|[i/at i-i/jv] = — 2(r+-f), 
6(0) = -2r, 6(1) = -2 7 , i.e., 



/ -2r 
2t 




(7) x(t) 



2 7 
-2(t + 7 ) 2 7 
2t 6 2 



V 







-2(7- - 

2t 


•<JV+1 



7) 2 7 
-2(t + 7 ) 2 7 
2r -2 7 / 

Using (J2) the corresponding 



x(t) 



subject to the initial condition x(0) = v G 
approximating PDE is then given by: 

d t u(t, z) = (t + ~/)h 2 d zz u(t, z) + 2(7 — r)hd z u(t, z) 

fs; <j ftu(i, 0) = 2 7 ^u(t, 0) + 2(7 - r)«(i, 0) 

9 t u(t, 1) = -2rhd z u(t, 1) - 2(7 - r)u(t, 1). 

To illustrate the effectiveness of our method numerically, we consider the special 
case of r = 7 = a/2, leading to the simplified equations 

' dtu(t, z) — ah 2 d zz u(t, z) 
(9) I d t u(t,0) = ahd z u{t,0) 

k d t u(t, 1) = -uhd z u{t, 1), 

where the associated generator has all its eigenvalues in (—00, 0]. Wishing to apply 
the Fourier method, we look for the solution in the form 

00 

3=0 

It is enough to find the eigenfunctions e Xjt Wj(z). The PDE and the boundary 
conditions then yield the system of equations 

Aw = ah 2 w" 



Aio(0) = ttW(0) 
t Xw{l) = -ato'(l). 
The first equation yields 

Wj(z) = ci^\ j cos(ujjz/h) + c 2 ,Aj sin(cjjz/h), 



° 2 ji > 0. Substituting into the first boundary condition we obtain 



with Xj = 

— WjCi Aj = c 2,\j, allowing us to choose c\.\ i = 1 and hence write 

u)j(z) = cos(ttJjz/h) — uij smitOjZ /h). 
Now substituting into the second boundary condition we obtain 



tan 



2u)j 



1 



This has exactly one solution in each interval ((2j — l)/if , (2j + l)/if ) for j > 0. 
The constants Cj are determined by the initial condition 

00 

u(0,z) = J2cjWj(z). 
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Figure 1. The probability distribution xu(t), k = 0, 1,2, ...,N 
at time t = 100 obtained from system ([7]) (circles) and the solution 
z H> u(t, z) of the PDE at time t = 100 (continuous line), with 
initially 200 nodes in state 1 and with N = 1000, r = 0.5, 7 = 0.5. 

In Figure [T] the solution of system Q is compared to the solution of the PDE ([SJ 
when t = 7 = 0.5, the latter was plotted using the Fourier method with the first 40 
eigenfunctions. The first 40 eigenvalues were determined by using Newton's method 
within each interval given above, and then we solved equation ([TU]) restricted to the 
first 40 variables. We observed that on our desktop computer MATLAB needed 
15.719000 seconds to get the ODE solution at t = 100, while for the Fourier method 
0.016000 seconds were needed to solve the PDE. 

4.2. SIS disease transmission model on a complete graph. The second moti- 
vation of our study comes from epidemiology where a paradigm disease transmission 
model is the simple susceptible-infected-susceptible (SIS) model on a completely 
connected graph with N + 1 nodes, i.e. all individuals are connected to each other. 
From the disease dynamic viewpoint, each individual is either susceptible (S) or 
infected (J) - a susceptible one with k + 1 infected neighbours can be infected at 
rate (fcr) and the infected ones can recover at a given rate (7) and become suscep- 
tible again. Since the graph is complete, the state space is the set {0, 1,2,..., N}, 
where a number k represents the state in which there are k infected nodes. Starting 
from state k the system can move either to state k + 1 or to k — 1 , since at a given 
instant only one node can change its state. When the system moves from state k to 
k + 1 then a susceptible node becomes infected. Hence the rate of this transition is 
k(N — fc)r, expressing that any of the N — k susceptible nodes can become infected 
and each of them has k infected neighbours (since the graph is complete). The 
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rate of transition from state k to k — 1 is Icy, because any of the k infected nodes 
can recover. Let us denote by Xk{t) the probability that the system is in state fc, 
i.e. there are k infected nodes. The above transition rates lead to the differential 
equation 



x{t) = (fc - 1)(N -k+ l)rx fc _i(t) - (k(N - k)r + ki)x{i) + (k + l)jx k+1 (t). 



(For k — and for k — N the equations contain only two terms.) Thus our system of 
ODEs takes the form given in (TI| with a k = k(N — k)r, c k — &7 and b k = —a k — c k , 
that is a(z) — N 2 tz{1 — z), c(z) — N^z and b(z) = —a(z) — c(z). We note that an 
approximation of this system by a first order PDE was investigated in Batkai, Kiss, 
Sikolya and Simon [5]. According to © our method yields the following second 
order approximation 



Our theorem implies that the solution of this PDE approximates the solution of 
the corresponding ODE ([1} in the order of 1/N 2 . We note that the usually used 
first order PDE approximates the ODE in the order of 1/N. The advantage of 
that first order PDE is that it can be solved analytically yielding the well-known 
mean-field approximation for the expected number of infected nodes. Our second 
order PDE cannot be solved analytically, hence only a numerical approximation can 
be obtained by using our method. It is the subject of future work to derive PDE 
approximations for epidemic propagation on different random graphs and compare 
their solutions to those of the original ODE system. 
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